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Introduction 


This  investigation  was  prompted  by  a  previous  report1,  describing 
airblast  calculations  made  with  the  HULL  hydrocode  for  the  MINERAL  FIND  2 
explosion  test,  involving  a  sphere  of  radius  38  cm,  packed  with  the 
explosive  HMX  and  detonated  at  304.8  cm  above  ground  (Figure  1).  A  plot 
of  pressure  contours  and  velocity  vectors  (Figure  2)  show  a  peculiar 
jetlike  behavior  in  pressure  and  velocity  at  t  -  0.75  msec.  The  HULL 
hydrocode  is  a  versatile  and  robust  code,  widely  used  for  compressible 
flow  calculations.  Based  on  previous  experience,  the  jet- like  behavior 
is  certainly  not  expected.  At  first  glance,  boundary  conditions  are 
suspected  to  be  the  primary  cause  of  this  anomaly. 

This  report  investigates  the  governing  equations  and  boundary 
conditions  that  are  used  in  the  HYDRO  subroutine  in  the  Eulerian, 
two-dimensional  calculation  mode  of  HULL.  It  was  our  primary  focus  given 
the  limited  time  of  the  investigation.  In  this  respect,  velocity, 
pressure  and  stress  boundary  conditions  were  checked  and,  if  necessary, 
modified.  It  was  anticipated  that  correcting  errors  in  the  boundary 
conditions  would  improve  the  time  and  accuracy  of  calculations. 

The  investigation  used  a  constant  subgrid  of  mesh  width  0.5  cm 
between  x  -  0  and  35  cm  for  all  trial  runs.  This  overcomes  the  grid 
dependency  of  the  scheme  developed.  First  velocity,  then  stress  boundary 
conditions  were  modified  and  tested.  Eventually  the  governing  equations 
had  to  be  reprogrammed  to  produce  successful  results.  The  following 
sections  present  the  analytical  theory,  a  discussion  of  the  computational 
and  experimental  results,  and  concluding  remarks.  The  appendix  shows  the 
SAIL  commands  implemented  in  this  investigation. 


Analytical  Formulation 


The  analytical  formulation  of  the  HULL  hydrocode  is  based  on  a 
mixture  of  Particle  in  Cell  (PIC)  and  Fluid  in  Cell  (FLIC)  codes 
originally  developed  at  the  Los  Alamos  Scientific  Laboratory  in  New 
Mexico  by  the  T-3  group.  The  code  is  very  large.  It  is  unnecessary  to 
report  all  the  details  here.  We  shall,  however,  focus  on  the  governing 
equations  and  boundary  conditions  in  the  subroutine  HYDRO.  For 
axisymmetric  flows,  the  continuity  and  momentum  equations  are  given  by 


1  C.  W.  Mastin,  "Mineral  Find  2  Airblast  Calculations",  Report 
presented  at  the  Airblast  Calculations  Meeting,  U.S.  Army  Engineer 
Waterways  Experiment  Station,  Vicksburg,  MS,  16  January,  1991, 
pp.  B  1-18. 
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where,  Vr  and  Vc  are  velocity  components  in  radial  and  axial  directions, 
Br  and  Bc  represent  body  force  per  unit  mass  along  r  and  z,  p  is  the 
density  and  D/Dt  has  the  usual  meaning  of  the  substantial  derivative. 

The  energy  equation  can  be  written  similarly  using  velocity  and  stress 
components . 

In  the  above  set  of  equations,  (2)  and  (3)  are  called  the  Navier's 
equations.  They  are  usually  cast  into  the  Navier-Stokes  form  using  the 
Stokes'  hypothesis  relating  stresses  to  the  strain  rates,  before  the 
above  set  can  be  solved  for  velocity  components.  In  the  HULL  context, 
however,  we  shall  retain  them  in  their  stress  forms. 

Whether  we  solve  the  governing  equations  in  the  Navier-Stokes  form 
or  retain  Navier's  form,  we  have  to  specify  boundary  conditions  in 
velocity  and  stress  components  along  the  boundaries  of  the  physical 
domain,  and,  as  in  the  context  of  a  time  dependent  problem,  specify 
initial  conditions.  Of  particular  interest  are  the  boundary  conditions 
along  the  line  of  symmetry  in  axisymmetric  calculations. 

There  are  several  possible  ways  the  governing  equations  in  the 
finite  difference  form.  The  HULL  hydrocode  utilizes  the  explicit 
calculation  mode,  where  velocity  components  are  expressed  at  the  edges  ox 
a  computational  cell  and  pressure  and  stress  components  are  expressed  at 
the  cell  center. 

Though  explicit  difference  schemes  have  the  drawback  of  small  time 
steps,  this  particular  code  has  superior  stability  properties  due  to  the 
"zip"  type  differencing  utilized  in  the  finite  differencing  mode.  This 
feature  also  maintains  the  second  order  accuracy,  even  with  linear 
interpolations.  The  details  of  the  donor  cell  differencing  are  widely 
available  in  literature  and  therefore  are  not  presented  here. 

The  calculation  in  HULL  proceeds  in  two  steps,  with  intermediate 
values  necessary  at  time  step  (n+h) .  In  the  first  phase  of  calculation, 
convection  in  terms  of  particle  transport  is  not  allowed.  The  finite 
difference  expressions  are  written  for  the  density  and  velocity  (r  and  z 
components)  equations,  (1)  through  (3).  In  addition,  the  finite 
difference  forms  for  the  specific  energy  and  pressure  equations  are 
derived  from: 
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In  Che  above  expressions,  e  is  the  specific  energy,  c  is  the  local 
speed  of  sound  in  the  medi.iim  and  p  is  the  pressure.  The  other  variables 
have  the  usual  meanings  as  before.  In  addition  to  the  above  equations, 
suitable  equations  of  state  are  used  depending  on  the  problem. 

At  this  stage,  equations  (1)  through  (5)  must  be  cast  into  finite 
difference  forms.  Although  explicit  differencing  can  be  used,  self 
adjoint  or  conservative  forms  must  be  maintained  in  differencing  for 
superior  error  properties.  This  is  the  reason  why  the  usual  flux 
differencing  of  radial  stress  terms  is  replaced  by  a  difference  between 
radial  and  hoop  stresses  in  quasi  one -dimensional  calculations.  The 
above  two  differencing  methods,  however,  do  not  differ  significantly  in 
axisymmetric  calculations  as  were  tested  here. 

The  above  equations  are  solved  with  the  help  of  two  types  of 
boundary  conditions  for  reflective  and  transmissive  boundaries .  In 
either  case,  additional  cell  values  for  fictitious  cells  are  calculated 
and  stored.  For  rigid  surfaces,  the  cells  adjacent  to  them  must  ensure 
no  flow  of  mass  or  energy  across  the  boundary.  Thus,  normal  velocity 
components  across  a  rigid  cell  must  be  zero  and  the  tangential  velocity 
components  must  be  preserved.  This  same  idea  pervades  in  normal  and 
tangential  stress  evaluations.  The  free  surfaces  must  be  free  of 
stresses  also.  Since  the  calculations  are  performed  for  inviscid  flows, 
the  boundary  conditions  at  a  rigid  surface  become  synonymous  with  those 
at  a  plane  of  symmetry.  To  keep  the  method  flexible,  similar 
differencing  as  in  the  equations  above,  using  hoop  stresses,  were 
developed  in  the  HULL  code. 


Numerical  Experiments 


A  thorough  visual  check  of  the  governing  equations  and  boundary 
conditions  was  made  prior  to  the  numerical  simulation.  The  size  of  HULL 
often  prohibits  detection  of  certain  anomalies  just  by  theoretical 
checks .  In  a  large  code  such  as  this ,  errors  do  happen  sometimes 
inadvertently.  However,  if  no  such  obvious  errors  are  found,  the  only 
means  to  identify  the  error  is  through  numerical  experiments.  One 
advantage  of  the  current  code  is  that  it  is  very  robust.  Minor  errors  by 
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the  pro  , rammer  are  often  forgiven.  However,  from  another  standpoint,  it 
beco  _;s  even  harder  to  make  the  program  fail  due  to  the  tweaking  of 
certain  inputs. 

In  performing  the  numerical  experiments ,  the  author  preferred  to 
use  a  theoretical  basis,  rather  than  rules  of  the  trade,  to  decide 
further  course  of  action.  For  example,  though  experience  shows  that  the 
use  of  the  difference  between  radial  and  hoop  stresses  may  be  substituted 
in  place  of  radial  conservative  differencing  (due  to  better  error 
cancellation),  this  is  not  the  general  practice  in  calculations  of  2-D 
cylindrical  geometry.  Therefore  an  attempt  was  made  to  first  cast  the 
formulation  the  way  it  should  be,  before  special  observations  such  as 
this  could  be  made. 

Some  features  of  the  HULL  calculation  scheme  were  retained  and  some 
modified  during  the  course  of  this  investigation.  For  example,  the  use 
of  the  minimum  value  of  pc2  was  retained  for  failure  calculations. 

However,  the  use  of  maximum  p  in  divergence  calculations  was  replaced  in 
some  runs  by  an  average  p. 

The  boundary  conditions  were  checked  first.  Apparently  there  were 
no  mistakes.  Several  cases  of  boundary  conditions  were  tested.  These 
included  intermediate  time  step  boundary  calculations  with  and  without 
particle  transport.  This  was  done  to  determine  the  sensitivity  of 
intermediate  boundary  values  on  the  final  results. 

In  some  experiments,  the  boundary  conditions  on  the  centerline  axis 
of  symmetry  were  cast  in  the  same  form  that  the  fundamental  differencing 
scheme  utilizes  to  treat  the  right  reflective  boundary.  It  should  be 
noted  that,  in  a  cell  calculation,  the  left  and  right  boundaries  of  the 
cell  must  conform  to  each  other  in  transporting  a  material  property. 

Finally,  when  no  significant  improvement  in  the  centerline  jet 
behavior  was  produced  by  the  use  of  modified  boundary  conditions, 
governing  equations  were  reprogrammed.  This  eventually  identified  the 
problem  in  the  specific  energy  equation  as  described  in  the  discussion  of 
results.  For  conciseness,  only  a  relevant  summary  will  be  presented  in 
the  next  section. 


Discussion  of  Results 


The  data  in  question  are  the  velocity  vector  and  pressure  plots 
along  the  centerline  in  axisymmetric  calculations.  The  original  data 
(reported  in  Reference  1  and  shown  in  Figure  2)  were  recreated  by  the 
author.  For  investigative  purposes,  all  test  runs  were  made  with  the 
airblast  calculations  of  the  same  HMX  detonation  problem.  This  first  set 
of  results  identifying  the  original  calculation  is  called  Test 
Calculation  No.  1,  and  shown  in  Figures  3  and  4.  Figure  3  shows  the 
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split  panel  plots  of  velocity  vectors  and  pressure  contours  at  time  t  - 
0.75  ms.  We  observe  spurious  jet  like  behavior  in  the  velocity  vector 
plot  along  the  centerline,  where  no  physical  influence  on  the  flow  is 
present.  The  pressure  contour  plots  on  the  left  panel  confirm  this 
behavior.  Figure  4  shows  the  corresponding  station  plot  at  ground  zero 
for  this  calculation.  Present  calculation  of  pressure  versus  time  is 
shown  by  the  solid  line.  This  plot  also  shows  two  broken  line  plots  of 
experimental  data  at  the  same  ground  zero  location.  The  shorter  broken 
lines  show  the  behavior  of  Airblast  Gage  No.  48  (AB-48)  data  (peak 
overpressure  of  93  MPa) ,  and  the  longer  broken  lines  show  the  behavior  of 
the  AB-47  data  (peak  overpressure  of  37  MPa).  For  comparison  purposes, 
all  test  calculations  will  be  presented  with  these  two  experimental  data 
plots.  The  wave  forms  generated  by  HULL  calculations  traditionally  lag 
the  experimental  data  in  arrival  time.  Therefore  the  abscissa  of  the 
calculational  results  has  been  shifted  by  200  /js  in  all  plots  of  this 
type. 


It  may  be  mentioned  that,  although  the  test  runs  were  made  from 
time  t  -  0  to  t  -  2.3  msec,  and  in  some  cases  up  to  t  -  3.3  msec,  for 
comparison  sake,  only  data  at  t  -  0.75  msec  will  be  compared.  In  the 
early  part  of  this  investigation,  temporal  growth  of  the  artifacts  were 
studied  at  different  times.  It  was  found  that  this  behavior  was  less 
prominent  after  the  shock  interacted  with  the  ground.  However,  in  some 
cases,  although  obscure,  the  centerline  vectors  show  up  in  a  single 
streak,  even  after  reflection  of  the  shock  from  the  ground. 

As  mentioned  before,  the  author  wished  to  determine  if  radial  flux 
differencing  would  improve  calculations  in  2-D  axially  symmetric 
problems .  With  flux  terms  cast  the  same  way  in  the  left  and  right 
boundaries,  the  results  produced  visible  improvements  in  velocity  vectors 
and  station  plots,  as  shown  in  Figures  5  and  6  (Test  Calculation  No.  2). 
This  input  uses  stress  values  on  the  left  boundary  set  at  zero.  The  peak 
overpressure  seems  to  attain  a  better  value  with  this  approach. 

The  best  peak  overpressure  was  obtained  in  the  next  case  (Figures  7 
and  8,  Test  Calculation  No.  3),  where  average  densities  are  used  instead 
of  maximum  densities  in  divergence  terms,  when  compared  with  the  previous 
test  case. 

Several  variations  of  governing  equation  and  boundary  conditions  in 
velocity  and  stresses  were  run.  However,  in  all  those  results,  velocity 
seemed  to  have  a  lingering  snag,  still  near  the  axis  of  symmetry  as 
apparent  in  both  the  Figures  5  and  7.  Therefore  those  additional  results 
are  not  presented  here. 

After  considerable  testing  failed  to  eliminate  the  jetlike  feature 
at  the  centerline,  focus  was  shifted  from  the  velocity  and  pressure  to 
the  specific  energy  calculations.  This  was  for  the  simple  reason  that 
although  pressure  looked  reasonable  (compare  Figure  3  with  Figures  5  or 
7),  velocity  still  seemed  to  be  having  an  adjustment  problem  from 
elsewhere,  which  points  to  the  specific  energy.  Finally,  it  turned  out 
that  only  those  changes  required  to  fix  the  artifacts  were  in  the 
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specific  energy  equation.  The  next  four  figures  (9  -  12)  show  the 
changes  implemented  in  the  radial  differencing  terms  of  the  specific 
energy  equation.  Figures  9  and  10  (Test  Calculation  No.  4)  are  for  this 
single  change  in  the  original  HULL  code.  Figures  11  and  12  (Test 
Calculation  No.  5)  reflect  the  case  of  radial  conservative  differencing 
in  velocity  equations  as  well.  Figures  9  and  11  both  show  correct  plots 
of  velocity  vectors  and  pressure.  However,  the  station  plots  have 
changed  in  these  runs  characteristically.  Now  the  peak  pressure  comes 
much  later.  At  this  time,  it  is  unknown  whether  the  fixing  of  the 
velocity  artifacts  has  more  merit  over  the  behavior  of  the  station  plots. 
This  phenomenon  needs  to  be  further  investigated  to  improve  HULL 
calculation  capabilities  for  investigating  this  type  of  explosion 
problem. 


Conclusions 


A  series  of  numerical  experiments  were  performed  to  determine  the 
cause  of  centerline  artifacts  in  HULL  axisymmetric  calculations. 

Although  it  appeared  originally  to  be  a  problem  of  velocity  or  pressure 
boundary  conditions,  the  jetlike  artifacts  were  found  to  be  caused  by  the 
specific  energy  equation.  At  this  stage,  the  centerline  instability 
problem  has  been  rectified.  However,  station  plots  show  that  earlier 
versions  appear  to  have  better  pressure  transients.  Further 
investigations  are  necessary  to  improve  the  pressure  calculations. 
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Figure  1.  Mineral  Find  2  test  configuration 
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Figure  2.  Original  Velocity  Vector  and  Pressure  Contour  Plot 
(reproduced  from  Reference  1) 
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Figure  4.  Station  Plot  at  Ground  Zero  for  Test  Calculation  No.  1 
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Figure  6.  Station  Plot  at  Ground  Zero  for  Test  Calculation  No.  2 
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Figure  7.  Split  Panel  Plots  for  Test  Calculation  No.  3 
(boundary  flux  plus  average  density  applied) 
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Figure  8.  Station  Plot  at  Ground  Zero  for  Test  Calculation  No.  3 
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Figure  9.  Split  Panel  Plots  for  Test  Calculation  No.  4 
(specific  energy  equation  modified) 
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Figure  10.  Station  Plot  at  Ground  Zero  for  Test  Calculation  No.  4 
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Figure  12.  Station  Plot  at  Ground  Zero  for  Test  Calculation  No.  5 
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Appendix:  SAIL  Input  for  Figures  3  through  12 

This  appendix  shows  the  sail  commands  for  the  run  of  each  of  the 
five  test  cases  presented  earlier  (Figures  3  through  12) .  Each  two 
consecutive  figures  were  obtained  as  a  result  of  a  set  of  new  SAIL 
commands  in  two  HULL  batch  run  and  starting  with  Figure  3.  Note  that 
SAIL  inputs  for  Test  Calculation  Nos.  2  through  5  show  only  the  changes 
from  the  SAIL  input  for  Test  Calculation  No.  1.  Other  trial  runs  that 
are  not  presented  here  may  be  obtained  from  the  author's  files  upon 
request. 
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c 

c 

c 

c********************  SAIL  input  for  Test  Calculation  Ho.  1  ************* 
c 
c 
c 

*d  50951 
cdir$  nolist 
*d  71495 

*skipto  *1  rezone 7  or  rezone9 
*i  86069 

$  disens ion  pres i(  i»ax_) ,presj (_jmax_) 

*i  86087 

data  omega,coef0,coefl,coef2,nosmoo/5. ,1. ,1. ,0. , 1/ 
c 

c  omega  -  ratio  of  maximum  to  minimum  grid  spacing 

c  coef0,l,2  -  coeficients  in  weight  function  for  adaptive  grid 

c  nosaoo  -  nuaber  of  smoothing  iterations  on  maxima  pressure 

c 

*i  86530 

*keepto  rez9  rezone9 
c 

c  this  is  an  adaptive  rezone  scheme  based  on  the  pressure, 

c  any  other  solution  variable  could  be  used  by  changing 

c  'n-l'  to  'n— variable  no  as  stored  on  tape4 ' . 

c 

j-0 

c 

c  read  in  pressure  and  compute  maxima  value  on  each  horizintal 

c  and  vertical  grid  line 

c 

do  19  ivy-1, nblks 
n-l 

do  11  jvy-l.nrawpb 
j-j+1 

*keepto  *2  not  incore 

nn-(jvy-l) *nvarpr 
call  bufin(j,h(nn+l) ,nvarpr) 
do  11  i-1, imax 
presi ( i ) -max (presi ( i) ,h(n)) 
pres j ( j ) -max (pres j ( j ), h (n) ) 
c 

c  use  these  instead  if  you  would  like  average  rather  than  maxima 

c  pressure 

c  presi (i) -presi (i)+h(n)*dy(j)/(y(jmax)-y(0)) 

c  presj (j)-presj (j)+h(n) *dx(i)/(x(imax)-x(0) ) 

c 

n-n+nh 

11  continue 

19  continue 

c 

c  smooth  the  one-dimensional  pressure  arrays  that  will  be  used 

c  to  define  the  weight  function 

c 

deltac-.25 

do  23  iterat-1 , nosaoo 
preso-presi(l) 
do  21  1-2, imax-1 
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prsave-pres i ( i ) 

deltax-2 . 0* (dx ( i+1) -dx ( i-1) ) / (dx ( i+1) +dx ( i-1) +2 . *dx ( i) ) 
deltat-min (0 . 5/  (abs  (del tax)  +1 .  e-4 ) , deltac) 
presi ( 1 ) —presi ( i)  + (presl (i+1) +preao-2 . +presi ( i) 

1  -. 5*deltax* (presi (i+1) -preso) ) *deltat 

preso-prsave 

21  continue 
preso-presj (1) 
do  22.  j«2 , jmax-1 
prsave— pres j ( j ) 

deltay-2.0*(dy(j+l)-dy(j-l))/(dy(j+l)+dy(j-l)+2.*dy(j)) 
deltat-min (0.5/ (abs (deltay) +1 . e-4 ) , deltac) 
presj (j)-presj (j)+(presj (j+l)+preso-2.*presj (j) 

1  5*deltay* (presj ( j+1) -preso) ) *deltat 

preso— prsave 

22  continue 

if (imax.gt.l)  then 
presi ( 1 ) -pres i ( 2 ) 
presi (imax) -presi ( imax-1) 
endif 

if (jmax.gt.l)  then 
presj (l)-presj (2) 
presj (j sax) -presj (j»ax-l) 
endif 

23  continue 
c 

c  smoothing  completed,  now  define  the  weight  function  and 

c  put  those  values  in  arrays  presi  and  presj 

c 

pmax-0 . 

preso-presi ( 1 ) 
do  12  i-2 , imax-1 
prsave— presi ( i) 

deltax-dx ( i+1) +2 . *dx ( i) +dx ( i-1 ) 

presi ( i) — coef 0+presi ( i) +coef l*abs (presi ( i+1) -preso) /deltax 

1  +4 . *coef 2*abs (presi ( i+1) +preso-2 . *presi ( i) 

2  - (dx( i+l) -dx (i-1) )* (presi (i+1) -preso) /deltax)/ (deltax) **2 
pmax— max (psax, presi (i) ) 

preso-prsave 

12  continue 

if (imax.gt. 1)  then 
presi ( imax) -presi ( imax-1) 
presi ( 1) -presi ( 2 ) 
endif 

preso-presj (1) 
do  13  j— 2, jmax-1 
prsave— pres j ( j ) 

deltay-dy ( j+1) +2 . *dy ( j ) +dy ( j-i) 

presj ( j ) -coef 0*pres j (j)+coefl*abs (presj (j+1) -preso) /deltay 

1  +4.*coef2*abs(presj ( j+1) +preso-2 . *presj (j) 

2  -(dy ( j+1) -dy ( j-l) ) * (presj (j+i) -preso) /deltay)/ (deltay) **2 
pmax-max(pmax, presj (j ) ) 

preso-prsave 

13  continue 

if (jmax.gt.l)  then 
presj (1) -presj (2) 
presj (jmax) -presj (jmax-1) 
endif 
c 

c  weight  function  is  modified  to  produce  desired  ratio  of 
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c  maxi  stun  to  minimum  grid  spacing 

c 

su*i«0. 

omemax— ( omega-1 . ) /pmax 

do  14  i-1, imax 

presi  ( i )  -l .  +omemax*presi  (  i ) 

14  sumi-sumi+i ./presi (i) 
sumj-O. 

do  15  j-l,jmax 

presj ( j ) — 1 . +omemax*pres j (j) 

15  sumj-sumj+1 . /presj ( j ) 
c 

c  new  grid  spacing  calculated;  satisfies:  dx*presi«constant 

c  and:  dy* presj -constant 

c 

xn(0)-x(0) 
yn(0)«y(0) 
do  16  i-1, imax 

dxn(i)-(x(imax)-x(0))/(presi(i) *sumi) 

xn(i)-xn(i-l)+dxn(i) 

rc (i) «xn ( i-1) + . 5*dxn (i) 

taun(i)-2.*pi*rc(i)*dxn(i) 

presi (i)-0. 

16  continue 

do  17  j-l,jmax 

«*yn  ( j )  -  (y  ( jMJt) -y  (0) )  /  (pres  j  ( j )  *su»j ) 

yn(j)*yn(j-l)+dyn(j) 

presj (j)-0. 

17  continue 
*label  rez9 

*i  120720 

common  /unit/nunit 
*i  121521 

nunit-66 

open  (nunit, file- 'plot. dat ' ) 

*i  130127 
c 

c  with  this  addition,  when  peel 1-. true. ,  the  grid  is  plotted 

c  (as  well  as  the  cell  indices  along  the  plot  border) . 

c  since  plotting  cell  indices  is  the  default,  you  must  use 

c  the  parameter  'nocells'  in  your  plot  statements  if  you  do 

c  not  want  a  grid:  for  example,  'pcont  (  nocells  )'. 

c  it  has  been  noted  that  when  pref-.true.,  the  grid  line  along 

c  the  axis  is  plotted  twice  and  that  grid  line  therefore  looks 

c  darker  than  the  other  lines, 

c 

if(pcell)  then 

ifirst-iminbd-1 

ilast-imaxbd 

j  f irst— jminbd-l 

jlast-jmaxbd 

xorg-0. 

if (pref .or. split)  xorg-xll 
if (split. and. left)  go  to  11 
do  9  i-if irst.  Hast 
xmesh-xorg+ ( xf ac*x ( i) -vxmin) /dvx 
ymesh- (yfac*y(j  first) -vymin) /dvy 
call  uplot(xmesh,ymesh,3) 
ymesh- (yfac*y (j last) -vymin) /dvy 
call  uplot ( xmesh , ymesh , 2 ) 
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9  continue 

do  10  j«j first, jlast 
ynesh- ( y f ac*y ( j ) -vymin) /dvy 
xMsh-xorg-f  (xf  ac*x  (  if  irst)  -vxmin)  /dvx 
call  upl ot  ( msb ,  ymesh ,  3 ) 

XM*h>xorg+  (xfac*x(ilast)  -vxain)  /dvx 
call  uplot(xmesh,yaesb,2) 

10  continue 

11  continue 

if (split. or.pref)  then 
if (split. and. .not. left)  go  to  14 
do  12  i-if irst,  Hast 
xmesh«Yorg- (xf ac*x ( i) -vxain) /dvx 
yaesh<-(yfac*y(  j  first)  -vy*in)/dvy 
call  uplot(xmesh,ymesh,3) 
ymesh— (yf ac*y ( j last) -vynin) /dvy 
call  uplot(xmesh,ymesh,2} 

12  continue 

do  13  j«j first, jlast 
ymesh-(yfac*y  ( j)  -vymin)/dvy 
xmesh—xorg- (xf ac*x ( if irst ) -vxmin) /dvx 
call  uplot (xmesh.ymesh, 3) 
xmesh— xorg- (xf ac*x ( ilast) -vxmin) /dvx 
call  uplot (xmesh, ymesh, 2) 

13  continue 

14  continue 
end  if 
endif 

*d  134919,134925 

*d  134925.2 
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SAIL  input  for  Test  Calculation  Mo.  2 


c 
c 
c 

C*********** 

c 
c 
c 

*d  5 0951 
cdir$  nolist 
*d  71495 
*skipto  *1  rezone7  or  rezone9 
*d  74708 

dv»2.*u(n) 

*d  74731 

srrl-0 . 
szzl-0. 
srzl-0. 

*d  74752 

ds»2 . *volf c*srz (n) *rc ( 1) /x ( 1 ) 

*d  74981 
*d  74983 

ds»(volfr*srr(nr) *rc(i+l) -volfe*srr(n) *rc(i) )/x(i) 

*d  75512 

1  (x(i)*srrr-x(i-l)*srrl) ) *vfac 
*i  86069 

$  distension  presi  (  imax  )  ,presj  (_jmax_) 

*i  86087 

data  oaega(coefO,coefl,coef2,nosaoo/5. ,1. ,1. ,0. ,1/ 
c 

c  omega  -  ratio  of  maximum  to  minimum  grid  spacing 

c  coef0,l,2  ■  coeficients  in  weight  function  for  adaptive  grid 

c  nosnoo  -  number  of  smoothing  iterations  on  maximum  pressure 

c 

*i  86530 

*keepto  rez9  rezone9 
c 

c  this  is  an  adaptive  rezone  scheme  based  on  the  pressure, 

c  any  other  solution  variable  could  be  used  by  changing 

c  'n-l'  to  'n-variable  no  as  stored  on  tape4 ' . 

c 

j-o 

c 

c  read  in  pressure  and  compute  maximum  value  on  each  horizintal 

c  and  vertical  grid  line 

c 


,etc. 
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c 

c 

c 

c********************  SAIL  input  for  Tost  calculation  Mo.  3  *********** 
c 
c 
c 

*d  50951 
cdir$  nolist 
*d  71495 

•skipto  *1  rezone7  or  rezone9 
*d  74416 

delp(c,d)«dt*dp/(0.5*(rhoc+rhoa)*(c+d-dt*dv) ) 

*d  74420 

dels(c,d)«dt*ds/  (0. 5*  (rhoc+rhoa)  *  (c+d-dt*dv) ) 

*d  74708 

dv»2.*u(n) 

*d  74731 

srrl-0. 

zzzl«0. 

*d  74752 

ds«2.*volfc*srz(n) *rc(l)/x(l) 

*d  74981 
*d  74983 

d*«* (volf r*srr (nr) *rc ( i+l ) -volf c*srr (n) *rc ( i ) )  /x  ( i) 

*d  75512 

1  (x(i) *srrr-x(i-l) *srrl) ) *vfac 

*i  86069 

$  dimension  presi(  imax  ) ,presj (_jsax_) 

*i  86087 

data  aaega,coefO,coefl,coef2,nosmoo/5.  ,1.  ,1.  ,0.  ,1/ 
c 

c  oMga  -  ratio  of  aaxiauH  to  minimus  grid  spacing 

c  coafO, 1,2  ■  cosficisnts  in  weight  function  for  adaptive  grid 

c  nosaoo  •  number  of  smooth ing  iterations  on  maximum  pressure 

c 

*i  86530 

*keepto  rez9  rezone9 
c 

c  this  is  an  adaptive  rezone  schene  based  on  the  pressure, 

c  any  other  solution  variable  could  be  used  by  changing 

c  'n-1'  to  'n-variable  no  as  stored  on  tape4 ' . 

c 

j-0 

c 

c  read  in  pressure  and  compute  maximum  value  on  each  horizintal 

c  and  vertical  grid  line 

c 


,etc. 
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c 

c 

c 

c********************  SAIL  input  for  Teat  Calculation  Ho.  4  **********’ 

c 

c 

c 

*d  50951 
cdir$  nolist 
*d  71495 

♦skipto  *l  rezone7  or  rezone9 
*d  75552 

xi(n)«xi(n)+dt*(2.0*pi*dy(j) *(ul*pl*x(i-l)-ur*pr*x(i) )/rc(i) 

*i  86069 

$  dimension  presi (imax  ) .presK  Imax  ) 

*i  86087 

data  omega,coef0,coefl,coef2,nosaoo/5.  ,1.  ,1.  ,0.  ,1/ 
c 

c  oaega  “  ratio  of  maxlaua  to  Minima  grid  spacing 

c  coafO, 1,2  -  coaficients  in  vaight  function  for  adaptive  grid 

c  nosaoo  -  number  of  smoothing  iterations  on  maximum  pressure 

c 

*i  86530 

*keepto  rez9  rezcne9 
c 

c  this  is  an  adaptive  rezone  scheme  based  on  the  pressure, 

c  any  other  solution  variable  could  be  used  by  changing 

c  'n-l'  to  #n«variable  no  as  stored  on  tape4 ' . 

c 

j-0 

c 

c  read  in  pressure  and  compute  maximum  value  on  each  horizintal 

c  and  vertical  grid  line 

c 


,etc. 
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c 

c 

c 

c******* *************  SAIL  input  for  Test  Calculation  No.  5  ********** 
c 
c 
c 

*d  S0951 
cdir$  nolist 

*d  714 95 

♦skip to  *l  rezone7  or  rezone9 
*d  74708 

dvn4.*u(n) 

*d  74731 

srrl-srr (n) 
szzl»o. 

*d  74752 

ds«4 . *volfc*srz (n) 

*d  74981 
*d  74983 

ds»  (volf  r*srr  (nr)  *rc  ( i+i)  -volf c*arr  (n)  *rc  (i))/x(i) 

*d  75512 

1  (x(i)*srrr-x(i-l)*srrl))*vfac 
*d  75552 

xi  (n)  «xi  (n)  +dt*  (2 . 0*pi*dy ( j )  *  (ul*pl*x  ( i-1)  -ur*pr*x( i) )  /rc(  i) 

*i  86069 

$  dimension  presi(  imax  ) .presj (_jmax_) 

*i  86087 

data  amega,coef0(coefl,coef2,nosmoo/5.,l.,l.  ,0.  ,1/ 
c 

c  oMga  -  ratio  of  maxima  to  minium  grid  spacing 

c  coef0,l,2  •  coeficients  in  weight  function  for  adaptive  grid 

c  nosmoo  •  number  of  smoothing  iterations  on  maximum  pressure 

c 

*i  86530 

♦keepto  rez9  rezone9 
c 

c  this  is  an  adaptive  rezone  scheme  based  on  the  pressure, 

c  any  other  solution  variable  could  be  used  by  changing 

c  'n-i'  to  '(invariable  no  as  stored  on  tape4 ' . 

c 

j-0 

c 

c  read  in  pressure  and  compute  maximum  value  on  each  horizintal 

c  and  vertical  grid  line 

c 


,etc 
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